home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
The Datafile PD-CD 1 Issue 2
/
PDCD-1 - Issue 02.iso
/
_utilities
/
utilities
/
001
/
meschach
/
!Meschach
/
c
/
torture
< prev
next >
Wrap
Text File
|
1994-08-04
|
28KB
|
1,058 lines
/**************************************************************************
**
** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
**
** Meschach Library
**
** This Meschach Library is provided "as is" without any express
** or implied warranty of any kind with respect to this software.
** In particular the authors shall not be liable for any direct,
** indirect, special, incidental or consequential damages arising
** in any way from use of the software.
**
** Everyone is granted permission to copy, modify and redistribute this
** Meschach Library, provided:
** 1. All copies contain this copyright notice.
** 2. All modified copies shall carry a notice stating who
** made the last modification and the date of such modification.
** 3. No charge is made for this software or works derived from it.
** This clause shall not be construed as constraining other software
** distributed on the same medium as this software, nor is a
** distribution fee considered a charge.
**
***************************************************************************/
/*
This file contains a series of tests for the Meschach matrix
library, parts 1 and 2
*/
static char rcsid[] = "$Id: torture.c,v 1.5 1994/05/18 01:53:58 des Exp $";
#include <stdio.h>
#include <math.h>
#include "matrix2.h"
#include "matlab.h"
#define errmesg(mesg) printf("Error: %s error: line %d\n",mesg,__LINE__)
#define notice(mesg) printf("# Testing %s...\n",mesg);
static char *test_err_list[] = {
"unknown error", /* 0 */
"testing error messages", /* 1 */
"unexpected end-of-file" /* 2 */
};
#define MAX_TEST_ERR (sizeof(test_err_list)/sizeof(char *))
/* extern int malloc_chain_check(); */
/* #define MEMCHK() if ( malloc_chain_check(0) ) \
{ printf("Error in malloc chain: \"%s\", line %d\n", \
__FILE__, __LINE__); exit(0); } */
#define MEMCHK()
/* cmp_perm -- returns 1 if pi1 == pi2, 0 otherwise */
int cmp_perm(pi1, pi2)
PERM *pi1, *pi2;
{
int i;
if ( ! pi1 || ! pi2 )
error(E_NULL,"cmp_perm");
if ( pi1->size != pi2->size )
return 0;
for ( i = 0; i < pi1->size; i++ )
if ( pi1->pe[i] != pi2->pe[i] )
return 0;
return 1;
}
/* px_rand -- generates sort-of random permutation */
PERM *px_rand(pi)
PERM *pi;
{
int i, j, k;
if ( ! pi )
error(E_NULL,"px_rand");
for ( i = 0; i < 3*pi->size; i++ )
{
j = (rand() >> 8) % pi->size;
k = (rand() >> 8) % pi->size;
px_transp(pi,j,k);
}
return pi;
}
#ifdef RISC_OS
#define SAVE_FILE "asx521_mat"
#else
#define SAVE_FILE "asx5213a.mat"
#endif
#define MATLAB_NAME "alpha"
char name[81] = MATLAB_NAME;
int main(argc, argv)
int argc;
char *argv[];
{
VEC *x = VNULL, *y = VNULL, *z = VNULL, *u = VNULL, *v = VNULL,
*w = VNULL;
VEC *diag = VNULL, *beta = VNULL;
PERM *pi1 = PNULL, *pi2 = PNULL, *pi3 = PNULL, *pivot = PNULL,
*blocks = PNULL;
MAT *A = MNULL, *Atmp = MNULL, *B = MNULL, *C = MNULL, *D = MNULL,
*Q = MNULL, *Qtmp = MNULL, *U = MNULL;
BAND *bA, *bB, *bC;
Real cond_est, s1, s2, s3;
int i, j, seed;
FILE *fp;
char *cp;
mem_info_on(TRUE);
setbuf(stdout,(char *)NULL);
seed = 1111;
if ( argc > 2 )
{
printf("usage: %s [seed]\n",argv[0]);
exit(0);
}
else if ( argc == 2 )
sscanf(argv[1], "%d", &seed);
/* set seed for rand() */
smrand(seed);
mem_stat_mark(1);
/* print version information */
m_version();
printf("# grep \"^Error\" the output for a listing of errors\n");
printf("# Don't panic if you see \"Error\" appearing; \n");
printf("# Also check the reported size of error\n");
printf("# This program uses randomly generated problems and therefore\n");
printf("# may occasionally produce ill-conditioned problems\n");
printf("# Therefore check the size of the error compared with MACHEPS\n");
printf("# If the error is within 1000*MACHEPS then don't worry\n");
printf("# If you get an error of size 0.1 or larger there is \n");
printf("# probably a bug in the code or the compilation procedure\n\n");
printf("# seed = %d\n",seed);
printf("# Check: MACHEPS = %g\n",MACHEPS);
/* allocate, initialise, copy and resize operations */
/* VEC */
notice("vector initialise, copy & resize");
x = v_get(12);
y = v_get(15);
z = v_get(12);
v_rand(x);
v_rand(y);
z = v_copy(x,z);
if ( v_norm2(v_sub(x,z,z)) >= MACHEPS )
errmesg("VEC copy");
v_copy(x,y);
x = v_resize(x,10);
y = v_resize(y,10);
if ( v_norm2(v_sub(x,y,z)) >= MACHEPS )
errmesg("VEC copy/resize");
x = v_resize(x,15);
y = v_resize(y,15);
if ( v_norm2(v_sub(x,y,z)) >= MACHEPS )
errmesg("VEC resize");
/* MAT */
notice("matrix initialise, copy & resize");
A = m_get(8,5);
B = m_get(3,9);
C = m_get(8,5);
m_rand(A);
m_rand(B);
C = m_copy(A,C);
if ( m_norm_inf(m_sub(A,C,C)) >= MACHEPS )
errmesg("MAT copy");
m_copy(A,B);
A = m_resize(A,3,5);
B = m_resize(B,3,5);
if ( m_norm_inf(m_sub(A,B,C)) >= MACHEPS )
errmesg("MAT copy/resize");
A = m_resize(A,10,10);
B = m_resize(B,10,10);
if ( m_norm_inf(m_sub(A,B,C)) >= MACHEPS )
errmesg("MAT resize");
MEMCHK();
/* PERM */
notice("permutation initialise, inverting & permuting vectors");
pi1 = px_get(15);
pi2 = px_get(12);
px_rand(pi1);
v_rand(x);
px_vec(pi1,x,z);
y = v_resize(y,x->dim);
pxinv_vec(pi1,z,y);
if ( v_norm2(v_sub(x,y,z)) >= MACHEPS )
errmesg("PERMute vector");
pi2 = px_inv(pi1,pi2);
pi3 = px_mlt(pi1,pi2,PNULL);
for ( i = 0; i < pi3->size; i++ )
if ( pi3->pe[i] != i )
errmesg("PERM inverse/multiply");
/* testing catch() etc */
notice("error handling routines");
catch(E_NULL,
catchall(v_add(VNULL,VNULL,VNULL);
errmesg("tracecatch() failure"),
printf("# tracecatch() caught error\n");
error(E_NULL,"main"));
errmesg("catch() failure"),
printf("# catch() caught E_NULL error\n"));
/* testing attaching a new error list (error list 2) */
notice("attaching error lists");
printf("# IT IS NOT A REAL WARNING ... \n");
err_list_attach(2,MAX_TEST_ERR,test_err_list,TRUE);
if (!err_is_list_attached(2))
errmesg("attaching the error list 2");
ev_err(__FILE__,1,__LINE__,"main",2);
err_list_free(2);
if (err_is_list_attached(2))
errmesg("detaching the error list 2");
/* testing inner products and v_mltadd() etc */
notice("inner products and linear combinations");
u = v_get(x->dim);
v_rand(u);
v_rand(x);
v_resize(y,x->dim);
v_rand(y);
v_mltadd(y,x,-in_prod(x,y)/in_prod(x,x),z);
if ( fabs(in_prod(x,z)) >= MACHEPS*x->dim )
errmesg("v_mltadd()/in_prod()");
s1 = -in_prod(x,y)/(v_norm2(x)*v_norm2(x));
sv_mlt(s1,x,u);
v_add(y,u,u);
if ( v_norm2(v_sub(u,z,u)) >= MACHEPS*x->dim )
errmesg("sv_mlt()/v_norm2()");
#ifdef ANSI_C
v_linlist(u,x,s1,y,1.0,VNULL);
if ( v_norm2(v_sub(u,z,u)) >= MACHEPS*x->dim )
errmesg("v_linlist()");
#endif
#ifdef VARARGS
v_linlist(u,x,s1,y,1.0,VNULL);
if ( v_norm2(v_sub(u,z,u)) >= MACHEPS*x->dim )
errmesg("v_linlist()");
#endif
MEMCHK();
/* vector norms */
notice("vector norms");
x = v_resize(x,12);
v_rand(x);
for ( i = 0; i < x->dim; i++ )
if ( v_entry(x,i) >= 0.5 )
v_set_val(x,i,1.0);
else
v_set_val(x,i,-1.0);
s1 = v_norm1(x);
s2 = v_norm2(x);
s3 = v_norm_inf(x);
if ( fabs(s1 - x->dim) >= MACHEPS*x->dim ||
fabs(s2 - sqrt((Real)(x->dim))) >= MACHEPS*x->dim ||
fabs(s3 - 1.0) >= MACHEPS )
errmesg("v_norm1/2/_inf()");
/* test matrix multiply etc */
notice("matrix multiply and invert");
A = m_resize(A,10,10);
B = m_resize(B,10,10);
m_rand(A);
m_inverse(A,B);
m_mlt(A,B,C);
for ( i = 0; i < C->m; i++ )
m_set_val(C,i,i,m_entry(C,i,i)-1.0);
if ( m_norm_inf(C) >= MACHEPS*m_norm_inf(A)*m_norm_inf(B)*5 )
errmesg("m_inverse()/m_mlt()");
MEMCHK();
/* ... and transposes */
notice("transposes and transpose-multiplies");
m_transp(A,A); /* can do square matrices in situ */
mtrm_mlt(A,B,C);
for ( i = 0; i < C->m; i++ )
m_set_val(C,i,i,m_entry(C,i,i)-1.0);
if ( m_norm_inf(C) >= MACHEPS*m_norm_inf(A)*m_norm_inf(B)*5 )
errmesg("m_transp()/mtrm_mlt()");
m_transp(A,A);
m_transp(B,B);
mmtr_mlt(A,B,C);
for ( i = 0; i < C->m; i++ )
m_set_val(C,i,i,m_entry(C,i,i)-1.0);
if ( m_norm_inf(C) >= MACHEPS*m_norm_inf(A)*m_norm_inf(B)*5 )
errmesg("m_transp()/mmtr_mlt()");
sm_mlt(3.71,B,B);
mmtr_mlt(A,B,C);
for ( i = 0; i < C->m; i++ )
m_set_val(C,i,i,m_entry(C,i,i)-3.71);
if ( m_norm_inf(C) >= MACHEPS*m_norm_inf(A)*m_norm_inf(B)*5 )
errmesg("sm_mlt()/mmtr_mlt()");
m_transp(B,B);
sm_mlt(1.0/3.71,B,B);
MEMCHK();
/* ... and matrix-vector multiplies */
notice("matrix-vector multiplies");
x = v_resize(x,A->n);
y = v_resize(y,A->m);
z = v_resize(z,A->m);
u = v_resize(u,A->n);
v_rand(x);
v_rand(y);
mv_mlt(A,x,z);
s1 = in_prod(y,z);
vm_mlt(A,y,u);
s2 = in_prod(u,x);
if ( fabs(s1 - s2) >= (MACHEPS*x->dim)*x->dim )
errmesg("mv_mlt()/vm_mlt()");
mv_mlt(B,z,u);
if ( v_norm2(v_sub(u,x,u)) >= MACHEPS*m_norm_inf(A)*m_norm_inf(B)*5 )
errmesg("mv_mlt()/m_inverse()");
MEMCHK();
/* get/set row/col */
notice("getting and setting rows and cols");
x = v_resize(x,A->n);
y = v_resize(y,B->m);
x = get_row(A,3,x);
y = get_col(B,3,y);
if ( fabs(in_prod(x,y) - 1.0) >= MACHEPS*m_norm_inf(A)*m_norm_inf(B)*5 )
errmesg("get_row()/get_col()");
sv_mlt(-1.0,x,x);
sv_mlt(-1.0,y,y);
set_row(A,3,x);
set_col(B,3,y);
m_mlt(A,B,C);
for ( i = 0; i < C->m; i++ )
m_set_val(C,i,i,m_entry(C,i,i)-1.0);
if ( m_norm_inf(C) >= MACHEPS*m_norm_inf(A)*m_norm_inf(B)*5 )
errmesg("set_row()/set_col()");
MEMCHK();
/* matrix norms */
notice("matrix norms");
A = m_resize(A,11,15);
m_rand(A);
s1 = m_norm_inf(A);
B = m_transp(A,B);
s2 = m_norm1(B);
if ( fabs(s1 - s2) >= MACHEPS*A->m )
errmesg("m_norm1()/m_norm_inf()");
C = mtrm_mlt(A,A,C);
s1 = 0.0;
for ( i = 0; i < C->m && i < C->n; i++ )
s1 += m_entry(C,i,i);
if ( fabs(sqrt(s1) - m_norm_frob(A)) >= MACHEPS*A->m*A->n )
errmesg("m_norm_frob");
MEMCHK();
/* permuting rows and columns */
notice("permuting rows & cols");
A = m_resize(A,11,15);
B = m_resize(B,11,15);
pi1 = px_resize(pi1,A->m);
px_rand(pi1);
x = v_resize(x,A->n);
y = mv_mlt(A,x,y);
px_rows(pi1,A,B);
px_vec(pi1,y,z);
mv_mlt(B,x,u);
if ( v_norm2(v_sub(z,u,u)) >= MACHEPS*A->m )
errmesg("px_rows()");
pi1 = px_resize(pi1,A->n);
px_rand(pi1);
px_cols(pi1,A,B);
pxinv_vec(pi1,x,z);
mv_mlt(B,z,u);
if ( v_norm2(v_sub(y,u,u)) >= MACHEPS*A->n )
errmesg("px_cols()");
MEMCHK();
/* MATLAB save/load */
notice("MATLAB save/load");
A = m_resize(A,12,11);
if ( (fp=fopen(SAVE_FILE,"w")) == (FILE *)NULL )
printf("Cannot perform MATLAB save/load test\n");
else
{
m_rand(A);
m_save(fp, A, name);
fclose(fp);
if ( (fp=fopen(SAVE_FILE,"r")) == (FILE *)NULL )
printf("Cannot open save file \"%s\"\n",SAVE_FILE);
else
{
M_FREE(B);
B = m_load(fp,&cp);
if ( strcmp(name,cp) || m_norm1(m_sub(A,B,B)) >= MACHEPS*A->m )
errmesg("mload()/m_save()");
}
}
MEMCHK();
/* Now, onto matrix factorisations */
A = m_resize(A,10,10);
B = m_resize(B,A->m,A->n);
m_copy(A,B);
x = v_resize(x,A->n);
y = v_resize(y,A->m);
z = v_resize(z,A->n);
u = v_resize(u,A->m);
v_rand(x);
mv_mlt(B,x,y);
z = v_copy(x,z);
notice("LU factor/solve");
pivot = px_get(A->m);
LUfactor(A,pivot);
tracecatch(LUsolve(A,pivot,y,x),"main");
tracecatch(cond_est = LUcondest(A,pivot),"main");
printf("# cond(A) approx= %g\n", cond_est);
if ( v_norm2(v_sub(x,z,u)) >= MACHEPS*v_norm2(x)*cond_est)
{
errmesg("LUfactor()/LUsolve()");
printf("# LU solution error = %g [cf MACHEPS = %g]\n",
v_norm2(v_sub(x,z,u)), MACHEPS);
}
v_copy(y,x);
tracecatch(LUsolve(A,pivot,x,x),"main");
tracecatch(cond_est = LUcondest(A,pivot),"main");
if ( v_norm2(v_sub(x,z,u)) >= MACHEPS*v_norm2(x)*cond_est)
{
errmesg("LUfactor()/LUsolve()");
printf("# LU solution error = %g [cf MACHEPS = %g]\n",
v_norm2(v_sub(x,z,u)), MACHEPS);
}
vm_mlt(B,z,y);
v_copy(y,x);
tracecatch(LUTsolve(A,pivot,x,x),"main");
if ( v_norm2(v_sub(x,z,u)) >= MACHEPS*v_norm2(x)*cond_est)
{
errmesg("LUfactor()/LUTsolve()");
printf("# LU solution error = %g [cf MACHEPS = %g]\n",
v_norm2(v_sub(x,z,u)), MACHEPS);
}
MEMCHK();
/* QR factorisation */
m_copy(B,A);
mv_mlt(B,z,y);
notice("QR factor/solve:");
diag = v_get(A->m);
beta = v_get(A->m);
QRfactor(A,diag);
QRsolve(A,diag,y,x);
if ( v_norm2(v_sub(x,z,u)) >= MACHEPS*v_norm2(x)*cond_est )
{
errmesg("QRfactor()/QRsolve()");
printf("# QR solution error = %g [cf MACHEPS = %g]\n",
v_norm2(v_sub(x,z,u)), MACHEPS);
}
Q = m_get(A->m,A->m);
makeQ(A,diag,Q);
makeR(A,A);
m_mlt(Q,A,C);
m_sub(B,C,C);
if ( m_norm1(C) >= MACHEPS*m_norm1(Q)*m_norm1(B) )
{
errmesg("QRfactor()/makeQ()/makeR()");
printf("# QR reconstruction error = %g [cf MACHEPS = %g]\n",
m_norm1(C), MACHEPS);
}
MEMCHK();
/* now try with a non-square matrix */
A = m_resize(A,15,7);
m_rand(A);
B = m_copy(A,B);
diag = v_resize(diag,A->n);
beta = v_resize(beta,A->n);
x = v_resize(x,A->n);
y = v_resize(y,A->m);
v_rand(y);
QRfactor(A,diag);
x = QRsolve(A,diag,y,x);
/* z is the residual vector */
mv_mlt(B,x,z); v_sub(z,y,z);
/* check B^T.z = 0 */
vm_mlt(B,z,u);
if ( v_norm2(u) >= MACHEPS*m_norm1(B)*v_norm2(y) )
{
errmesg("QRfactor()/QRsolve()");
printf("# QR solution error = %g [cf MACHEPS = %g]\n",
v_norm2(u), MACHEPS);
}
Q = m_resize(Q,A->m,A->m);
makeQ(A,diag,Q);
makeR(A,A);
m_mlt(Q,A,C);
m_sub(B,C,C);
if ( m_norm1(C) >= MACHEPS*m_norm1(Q)*m_norm1(B) )
{
errmesg("QRfactor()/makeQ()/makeR()");
printf("# QR reconstruction error = %g [cf MACHEPS = %g]\n",
m_norm1(C), MACHEPS);
}
D = m_get(A->m,Q->m);
mtrm_mlt(Q,Q,D);
for ( i = 0; i < D->m; i++ )
m_set_val(D,i,i,m_entry(D,i,i)-1.0);
if ( m_norm1(D) >= MACHEPS*m_norm1(Q)*m_norm_inf(Q) )
{
errmesg("QRfactor()/makeQ()/makeR()");
printf("# QR orthogonality error = %g [cf MACHEPS = %g]\n",
m_norm1(D), MACHEPS);
}
MEMCHK();
/* QRCP factorisation */
m_copy(B,A);
notice("QR factor/solve with column pivoting");
pivot = px_resize(pivot,A->n);
QRCPfactor(A,diag,pivot);
z = v_resize(z,A->n);
QRCPsolve(A,diag,pivot,y,z);
/* pxinv_vec(pivot,z,x); */
/* now compute residual (z) vector */
mv_mlt(B,x,z); v_sub(z,y,z);
/* check B^T.z = 0 */
vm_mlt(B,z,u);
if ( v_norm2(u) >= MACHEPS*m_norm1(B)*v_norm2(y) )
{
errmesg("QRCPfactor()/QRsolve()");
printf("# QR solution error = %g [cf MACHEPS = %g]\n",
v_norm2(u), MACHEPS);
}
Q = m_resize(Q,A->m,A->m);
makeQ(A,diag,Q);
makeR(A,A);
m_mlt(Q,A,C);
M_FREE(D);
D = m_get(B->m,B->n);
px_cols(pivot,C,D);
m_sub(B,D,D);
if ( m_norm1(D) >= MACHEPS*m_norm1(Q)*m_norm1(B) )
{
errmesg("QRCPfactor()/makeQ()/makeR()");
printf("# QR reconstruction error = %g [cf MACHEPS = %g]\n",
m_norm1(D), MACHEPS);
}
MEMCHK();
/* Cholesky and LDL^T factorisation */
/* Use these for normal equations approach */
notice("Cholesky factor/solve");
mtrm_mlt(B,B,A);
CHfactor(A);
u = v_resize(u,B->n);
vm_mlt(B,y,u);
z = v_resize(z,B->n);
CHsolve(A,u,z);
v_sub(x,z,z);
if ( v_norm2(z) >= MACHEPS*v_norm2(x)*100 )
{
errmesg("CHfactor()/CHsolve()");
printf("# Cholesky solution error = %g [cf MACHEPS = %g]\n",
v_norm2(z), MACHEPS);
}
/* modified Cholesky factorisation should be identical with Cholesky
factorisation provided the matrix is "sufficiently positive definite */
mtrm_mlt(B,B,C);
MCHfactor(C,MACHEPS);
m_sub(A,C,C);
if ( m_norm1(C) >= MACHEPS*m_norm1(A) )
{
errmesg("MCHfactor()");
printf("# Modified Cholesky error = %g [cf MACHEPS = %g]\n",
m_norm1(C), MACHEPS);
}
/* now test the LDL^T factorisation -- using a negative def. matrix */
mtrm_mlt(B,B,A);
sm_mlt(-1.0,A,A);
m_copy(A,C);
LDLfactor(A);
LDLsolve(A,u,z);
w = v_get(A->m);
mv_mlt(C,z,w);
v_sub(w,u,w);
if ( v_norm2(w) >= MACHEPS*v_norm2(u)*m_norm1(C) )
{
errmesg("LDLfactor()/LDLsolve()");
printf("# LDL^T residual = %g [cf MACHEPS = %g]\n",
v_norm2(w), MACHEPS);
}
v_add(x,z,z);
if ( v_norm2(z) >= MACHEPS*v_norm2(x)*100 )
{
errmesg("LDLfactor()/LDLsolve()");
printf("# LDL^T solution error = %g [cf MACHEPS = %g]\n",
v_norm2(z), MACHEPS);
}
MEMCHK();
/* and now the Bunch-Kaufman-Parlett method */
/* set up D to be an indefinite diagonal matrix */
notice("Bunch-Kaufman-Parlett factor/solve");
D = m_resize(D,B->m,B->m);
m_zero(D);
w = v_resize(w,B->m);
v_rand(w);
for ( i = 0; i < w->dim; i++ )
if ( v_entry(w,i) >= 0.5 )
m_set_val(D,i,i,1.0);
else
m_set_val(D,i,i,-1.0);
/* set A <- B^T.D.B */
C = m_resize(C,B->n,B->n);
C = mtrm_mlt(B,D,C);
A = m_mlt(C,B,A);
C = m_resize(C,B->n,B->n);
C = m_copy(A,C);
/* ... and use BKPfactor() */
blocks = px_get(A->m);
pivot = px_resize(pivot,A->m);
x = v_resize(x,A->m);
y = v_resize(y,A->m);
z = v_resize(z,A->m);
v_rand(x);
mv_mlt(A,x,y);
BKPfactor(A,pivot,blocks);
printf("# BKP pivot =\n"); px_output(pivot);
printf("# BKP blocks =\n"); px_output(blocks);
BKPsolve(A,pivot,blocks,y,z);
/* compute & check residual */
mv_mlt(C,z,w);
v_sub(w,y,w);
if ( v_norm2(w) >= MACHEPS*m_norm1(C)*v_norm2(z) )
{
errmesg("BKPfactor()/BKPsolve()");
printf("# BKP residual size = %g [cf MACHEPS = %g]\n",
v_norm2(w), MACHEPS);
}
/* check update routines */
/* check LDLupdate() first */
notice("update L.D.L^T routine");
A = mtrm_mlt(B,B,A);
m_resize(C,A->m,A->n);
C = m_copy(A,C);
LDLfactor(A);
s1 = 3.7;
w = v_resize(w,A->m);
v_rand(w);
for ( i = 0; i < C->m; i++ )
for ( j = 0; j < C->n; j++ )
m_set_val(C,i,j,m_entry(C,i,j)+s1*v_entry(w,i)*v_entry(w,j));
LDLfactor(C);
LDLupdate(A,w,s1);
/* zero out strictly upper triangular parts of A and C */
for ( i = 0; i < A->m; i++ )
for ( j = i+1; j < A->n; j++ )
{
m_set_val(A,i,j,0.0);
m_set_val(C,i,j,0.0);
}
if ( m_norm1(m_sub(A,C,C)) >= sqrt(MACHEPS)*m_norm1(A) )
{
errmesg("LDLupdate()");
printf("# LDL update matrix error = %g [cf MACHEPS = %g]\n",
m_norm1(C), MACHEPS);
}
/* BAND MATRICES */
#define COL 40
#define UDIAG 5
#define LDIAG 2
smrand(101);
bA = bd_get(LDIAG,UDIAG,COL);
bB = bd_get(LDIAG,UDIAG,COL);
bC = bd_get(LDIAG,UDIAG,COL);
A = m_resize(A,COL,COL);
B = m_resize(B,COL,COL);
pivot = px_resize(pivot,COL);
x = v_resize(x,COL);
w = v_resize(w,COL);
z = v_resize(z,COL);
m_rand(A);
/* generate band matrix */
mat2band(A,LDIAG,UDIAG,bA);
band2mat(bA,A); /* now A is banded */
bB = bd_copy(bA,bB);
v_rand(x);
mv_mlt(A,x,w);
z = v_copy(w,z);
notice("band LU factorization");
bdLUfactor(bA,pivot);
/* pivot will be changed */
bdLUsolve(bA,pivot,z,z);
v_sub(x,z,z);
if (v_norm2(z) > v_norm2(x)*sqrt(MACHEPS)) {
errmesg("incorrect solution (band LU factorization)");
printf(" ||exact sol. - computed sol.|| = %g [MACHEPS = %g]\n",
v_norm2(z),MACHEPS);
}
/* solve transpose system */
notice("band LU factorization for transpose system");
m_transp(A,B);
mv_mlt(B,x,w);
bd_copy(bB,bA);
bd_transp(bA,bA);
/* transposition in situ */
bd_transp(bA,bA);
bd_transp(bA,bB);
bdLUfactor(bB,pivot);
bdLUsolve(bB,pivot,w,z);
v_sub(x,z,z);
if (v_norm2(z) > v_norm2(x)*sqrt(MACHEPS)) {
errmesg("incorrect solution (band transposed LU factorization)");
printf(" ||exact sol. - computed sol.|| = %g [MACHEPS = %g]\n",
v_norm2(z),MACHEPS);
}
/* Cholesky factorization */
notice("band Choleski LDL' factorization");
m_add(A,B,A); /* symmetric matrix */
for (i=0; i < COL; i++) /* positive definite */
A->me[i][i] += 2*LDIAG;
mat2band(A,LDIAG,LDIAG,bA);
band2mat(bA,A); /* corresponding matrix A */
v_rand(x);
mv_mlt(A,x,w);
z = v_copy(w,z);
bdLDLfactor(bA);
z = bdLDLsolve(bA,z,z);
v_sub(x,z,z);
if (v_norm2(z) > v_norm2(x)*sqrt(MACHEPS)) {
errmesg("incorrect solution (band LDL' factorization)");
printf(" ||exact sol. - computed sol.|| = %g [MACHEPS = %g]\n",
v_norm2(z),MACHEPS);
}
/* new bandwidths */
m_rand(A);
bA = bd_resize(bA,UDIAG,LDIAG,COL);
bB = bd_resize(bB,UDIAG,LDIAG,COL);
mat2band(A,UDIAG,LDIAG,bA);
band2mat(bA,A);
bd_copy(bA,bB);
mv_mlt(A,x,w);
notice("band LU factorization (resized)");
bdLUfactor(bA,pivot);
/* pivot will be changed */
bdLUsolve(bA,pivot,w,z);
v_sub(x,z,z);
if (v_norm2(z) > v_norm2(x)*sqrt(MACHEPS)) {
errmesg("incorrect solution (band LU factorization)");
printf(" ||exact sol. - computed sol.|| = %g [MACHEPS = %g]\n",
v_norm2(z),MACHEPS);
}
/* testing transposition */
notice("band matrix transposition");
m_zero(bA->mat);
bd_copy(bB,bA);
m_zero(bB->mat);
bd_copy(bA,bB);
bd_transp(bB,bB);
bd_transp(bB,bB);
m_zero(bC->mat);
bd_copy(bB,bC);
m_sub(bA->mat,bC->mat,bC->mat);
if (m_norm_inf(bC->mat) > MACHEPS*bC->mat->n) {
errmesg("band transposition");
printf(" difference ||A - (A')'|| = %g\n",m_norm_inf(bC->mat));
}
bd_free(bA);
bd_free(bB);
bd_free(bC);
MEMCHK();
/* now check QRupdate() routine */
notice("update QR routine");
B = m_resize(B,15,7);
A = m_resize(A,B->m,B->n);
m_copy(B,A);
diag = v_resize(diag,A->n);
beta = v_resize(beta,A->n);
QRfactor(A,diag);
Q = m_resize(Q,A->m,A->m);
makeQ(A,diag,Q);
makeR(A,A);
m_resize(C,A->m,A->n);
w = v_resize(w,A->m);
v = v_resize(v,A->n);
u = v_resize(u,A->m);
v_rand(w);
v_rand(v);
vm_mlt(Q,w,u);
QRupdate(Q,A,u,v);
m_mlt(Q,A,C);
for ( i = 0; i < B->m; i++ )
for ( j = 0; j < B->n; j++ )
m_set_val(B,i,j,m_entry(B,i,j)+v_entry(w,i)*v_entry(v,j));
m_sub(B,C,C);
if ( m_norm1(C) >= MACHEPS*m_norm1(A)*m_norm1(Q)*2 )
{
errmesg("QRupdate()");
printf("# Reconstruction error in QR update = %g [cf MACHEPS = %g]\n",
m_norm1(C), MACHEPS);
}
m_resize(D,Q->m,Q->n);
mtrm_mlt(Q,Q,D);
for ( i = 0; i < D->m; i++ )
m_set_val(D,i,i,m_entry(D,i,i)-1.0);
if ( m_norm1(D) >= 10*MACHEPS*m_norm1(Q)*m_norm_inf(Q) )
{
errmesg("QRupdate()");
printf("# QR update orthogonality error = %g [cf MACHEPS = %g]\n",
m_norm1(D), MACHEPS);
}
/* Now check eigenvalue/SVD routines */
notice("eigenvalue and SVD routines");
A = m_resize(A,11,11);
B = m_resize(B,A->m,A->n);
C = m_resize(C,A->m,A->n);
D = m_resize(D,A->m,A->n);
Q = m_resize(Q,A->m,A->n);
m_rand(A);
/* A <- A + A^T for symmetric case */
m_add(A,m_transp(A,C),A);
u = v_resize(u,A->m);
u = symmeig(A,Q,u);
m_zero(B);
for ( i = 0; i < B->m; i++ )
m_set_val(B,i,i,v_entry(u,i));
m_mlt(Q,B,C);
mmtr_mlt(C,Q,D);
m_sub(A,D,D);
if ( m_norm1(D) >= MACHEPS*m_norm1(Q)*m_norm_inf(Q)*v_norm_inf(u)*3 )
{
errmesg("symmeig()");
printf("# Reconstruction error = %g [cf MACHEPS = %g]\n",
m_norm1(D), MACHEPS);
}
mtrm_mlt(Q,Q,D);
for ( i = 0; i < D->m; i++ )
m_set_val(D,i,i,m_entry(D,i,i)-1.0);
if ( m_norm1(D) >= MACHEPS*m_norm1(Q)*m_norm_inf(Q)*3 )
{
errmesg("symmeig()");
printf("# symmeig() orthogonality error = %g [cf MACHEPS = %g]\n",
m_norm1(D), MACHEPS);
}
MEMCHK();
/* now test (real) Schur decomposition */
/* m_copy(A,B); */
M_FREE(A);
A = m_get(11,11);
m_rand(A);
B = m_copy(A,B);
MEMCHK();
B = schur(B,Q);
MEMCHK();
m_mlt(Q,B,C);
mmtr_mlt(C,Q,D);
MEMCHK();
m_sub(A,D,D);
if ( m_norm1(D) >= MACHEPS*m_norm1(Q)*m_norm_inf(Q)*m_norm1(B)*5 )
{
errmesg("schur()");
printf("# Schur reconstruction error = %g [cf MACHEPS = %g]\n",
m_norm1(D), MACHEPS);
}
/* orthogonality check */
mmtr_mlt(Q,Q,D);
for ( i = 0; i < D->m; i++ )
m_set_val(D,i,i,m_entry(D,i,i)-1.0);
if ( m_norm1(D) >= MACHEPS*m_norm1(Q)*m_norm_inf(Q)*10 )
{
errmesg("schur()");
printf("# Schur orthogonality error = %g [cf MACHEPS = %g]\n",
m_norm1(D), MACHEPS);
}
MEMCHK();
/* check for termination */
Atmp = m_get(2,2);
Qtmp = m_get(2,2);
/* this is a 2 x 2 matrix with real eigenvalues */
Atmp->me[0][0] = 1;
Atmp->me[1][1] = 1;
Atmp->me[0][1] = 4;
Atmp->me[1][0] = 1;
schur(Atmp,Qtmp);
MEMCHK();
/* now test SVD */
A = m_resize(A,11,7);
m_rand(A);
U = m_get(A->n,A->n);
Q = m_resize(Q,A->m,A->m);
u = v_resize(u,max(A->m,A->n));
svd(A,Q,U,u);
/* check reconstruction of A */
D = m_resize(D,A->m,A->n);
C = m_resize(C,A->m,A->n);
m_zero(D);
for ( i = 0; i < min(A->m,A->n); i++ )
m_set_val(D,i,i,v_entry(u,i));
mtrm_mlt(Q,D,C);
m_mlt(C,U,D);
m_sub(A,D,D);
if ( m_norm1(D) >= MACHEPS*m_norm1(U)*m_norm_inf(Q)*m_norm1(A) )
{
errmesg("svd()");
printf("# SVD reconstruction error = %g [cf MACHEPS = %g]\n",
m_norm1(D), MACHEPS);
}
/* check orthogonality of Q and U */
D = m_resize(D,Q->n,Q->n);
mtrm_mlt(Q,Q,D);
for ( i = 0; i < D->m; i++ )
m_set_val(D,i,i,m_entry(D,i,i)-1.0);
if ( m_norm1(D) >= MACHEPS*m_norm1(Q)*m_norm_inf(Q)*5 )
{
errmesg("svd()");
printf("# SVD orthognality error (Q) = %g [cf MACHEPS = %g\n",
m_norm1(D), MACHEPS);
}
D = m_resize(D,U->n,U->n);
mtrm_mlt(U,U,D);
for ( i = 0; i < D->m; i++ )
m_set_val(D,i,i,m_entry(D,i,i)-1.0);
if ( m_norm1(D) >= MACHEPS*m_norm1(U)*m_norm_inf(U)*5 )
{
errmesg("svd()");
printf("# SVD orthognality error (U) = %g [cf MACHEPS = %g\n",
m_norm1(D), MACHEPS);
}
for ( i = 0; i < u->dim; i++ )
if ( v_entry(u,i) < 0 || (i < u->dim-1 &&
v_entry(u,i+1) > v_entry(u,i)) )
break;
if ( i < u->dim )
{
errmesg("svd()");
printf("# SVD sorting error\n");
}
/* test of long vectors */
notice("Long vectors");
x = v_resize(x,100000);
y = v_resize(y,100000);
z = v_resize(z,100000);
v_rand(x);
v_rand(y);
v_mltadd(x,y,3.0,z);
sv_mlt(1.0/3.0,z,z);
v_mltadd(z,x,-1.0/3.0,z);
v_sub(z,y,x);
if (v_norm2(x) >= MACHEPS*(x->dim)) {
errmesg("long vectors");
printf(" norm = %g\n",v_norm2(x));
}
mem_stat_free(1);
MEMCHK();
/**************************************************
VEC *x, *y, *z, *u, *v, *w;
VEC *diag, *beta;
PERM *pi1, *pi2, *pi3, *pivot, *blocks;
MAT *A, *B, *C, *D, *Q, *U;
**************************************************/
V_FREE(x); V_FREE(y); V_FREE(z);
V_FREE(u); V_FREE(v); V_FREE(w);
V_FREE(diag); V_FREE(beta);
PX_FREE(pi1); PX_FREE(pi2); PX_FREE(pi3);
PX_FREE(pivot); PX_FREE(blocks);
M_FREE(A); M_FREE(B); M_FREE(C);
M_FREE(D); M_FREE(Q); M_FREE(U);
M_FREE(Atmp); M_FREE(Qtmp);
MEMCHK();
printf("# Finished torture test\n");
mem_info();
return 0;
}